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Abstract 

In this note we present a flexible approach to perform the propagation of track 
parameters and their derivatives in an inhomogeneous magnetic field, keeping the 
computational effort small. We discuss also a Kalman filter implementation using 
this optimized computation of the derivatives. 

1 Introduction 

Dedicated hadronic B factory experiments as HERA-B are designed as forward spectrom- 
eters, in adjustment to the huge Lorentz boost under which the B decay particles are 
produced. The HERA-B tracking concept [1] relies on the propagation of track candi- 
dates which have been found in the pattern tracker in the field-free area upstream through 
the spectrometer magnet. The concurrent track evolution strategy [2] which employs the 
Kalman filter technique is used to cope with the large track densities and still give a 
high track finding efficiency. Practical implementation of this concept requires a fast 
and precise procedure to transport both the track parameters and its covariance matrix 
estimate, in spite of the inhomogeneity of the magnetic field. First solution based on a 
fifth-order Runge-Kutta method have been shown in [3]. 

This note presents the mathematical basis and the program implementation of an 
approach which was developed to achieve a further gain in speed, and at the same time 
warrant sufficient accuracy for an optimal operation of the track finding process. This 
was achieved by providing a set of methods which are optimized for different transport 
ranges, and by testing the validity of each approximation directly within the track finding 
application [4, 5]. 

The Kalman filter technique is used very often in high energy physics experiments. We 
include Sec. 4 in the revised version of the note to discuss the optimized implementation of 
the Kalman filter for magnet tracking, which reduces strongly the amount of computation. 
It was not described in detail in [4, 5]. 
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2 Equations of motion and solution 



In the following we will use a coordinate system in which the z axis points along the 
proton beam, the x axis is directed normal to it in the horizontal plane, pointing towards 
the inside of the HERA ring, and the y axis is oriented upwards such that x, y and z form 
a right-handed system. This system is identical to the Arte coordinate system defined in 
the appendix of [6]. 

The following choice of track parameters is suited for fixed target experiments with 
relatively small transverse momenta. A particle with momentum p, charge Q e and coor- 
dinates x is described by the state vector at a reference z coordinate: 



/ x \ 



x(z) 



by 
\ 1 J 



(1) 



where x and y are the transverse coordinates, t x = p x /p z and t y = p y /p z are the track 
slopes, and q = Q e / \p\. The parameter Q e is the particle charge in units of the elementary 
charge. We use the following units: x,y,z are in centimeters, p is in GeV/c and the 
magnetic field B in kGauss. 

The trajectory of a particle in a static magnetic field B(x), neglecting stochastic 
perturbations as energy loss and multiple scattering, must satisfy the equations of mo- 
tion [7] : 



dx/dz - 
dy/dz = 
dt x /dz 
dty/dz 
Q = 



t 



q ■ v ■ A x (t x , ty, 



B), (2) 

— * 

q-V ■ Ay(t X ,ty,B), 

where parameter v is proportional to the velocity of light and is therefore defined as 

0.000299792458 (GeV/c) kG' 1 



v 



and the functions A x ,A y are 
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The initial values at z = z are 



x o — ( x o, Do, t x o, t y o, qo). 



(3) 



Let us assume we should propagate the particle parameters from plane z to plane z. For 
solution of the equations (2) three different methods are used. The choice depends on 
the distance s = z — z Q between these planes. 

1. \s\ < 20 cm: a parabolic expansion of the particle trajectory is used 
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2. 20 cm < \s\ < 60 cm: the classical fourth-order Runge-Kutta method [8] is 
selected to find solution of the equations (2) . 

3. \s\ > 60 cm: a fifth-order Runge-Kutta method with adaptive step size control [8] 
is used. 

Use of the Kalman filter technique [9] for pattern recognition requires that particle 
parameters and their covariance matrix can be transported to the location of the next 
measurement. Evaluation of the derivatives of the state vector components with respect 
to their initial values (3) is needed to transport the covariance matrix. This is achieved 
by integrating the equations for derivatives together with the 'zero trajectory' (2). 



3 Equations for derivatives 

Let us assume the magnetic field is smooth enough and field gradients can be neglected. 
In this case, the equations (2) are invariant with respect to small shifts by x and y. This 
means that derivatives with respect to initial xo, yo are trivial : 

dx T /dx = (1, 0, 0, 0, 0) , 
dx T /dy = (0, 1, 0, 0, 0) . 

To obtain equations for dx/dt x o, let us differentiate equations (2) with respect to t x0 and 
change the order of the derivative operators d/dt x0 and d/dz on the left hand sides : 

d/dz(dx/dt xQ ) = dt x /dt x0 , 

d/dz(dy/dt x0 ) = dt y /dt x0 , 

d/dz(dt x /dt x0 ) = q -v [(dA x /dt x )(dt x /dt x0 ) + (dA x /dt y )(dt y /dt x0 )} , (4) 

d/dz(dt y /dt x0 ) = q -v ■ [(dA y /dt x )(dt x /dt x0 ) 4 (dA y /dt y )(dt y /dt x0 )] , 

dq/dt x0 = 0, 

where 

dA x /dt x = t x ■ A x /(1 + t 2 x + t 2 y ) + (1 + t 2 x + t 2 y y* ■ (t y ■ B x - 2 • t x ■ By) , 

dA X /dty = ty ■ A X /{\ +tl+ t 2 y ) + (1 + t\ + ■ (t x ■ B x + B z ) , 

dAy/8t X = t X ■ Ay /{I +tl + t\) + (1 + tl + t 2 y ) l 2 ■ (~ty • By ~ B g ) , 
dAy/dty = ty • A y / ( 1 + t\ 4" ^) 4" ( 1 4" t ^ + t J ) ^ • [~t X ' B 'y + 2 • ty • B X ) . 

Initial values for the solution of equations (4) are : 

dx T /dt x0 = (0, 0, 1, 0, 0). (5) 

The equations for dx/dt y o are similar to equations (4) , but the initial values are : 

dx T /dt y0 = (0, 0, 0, 1, 0) . 

To obtain equations for dx/dq , let us differentiate the equations (2) with respect to q 
and change the order of the derivative operators d/dqo and d/dz in the left parts : 

d/dz(dx/dq ) = dt x /dq , 

d/dz(dy/dq ) = dt y /dq , 

d/dz(dt x /dq ) = v-A x + vq - [(dA x /dt x )(dt x /dq ) 4 (dA x /dt y )(dt y /dq )} , (6) 

d/dz(dt y /dq ) = v ■ A y + v ■ q ■ [(dA y /dt x )(dt x /dq ) + (dA y /dt y )(dt y /dq )] , 

dq/dqo = 1 . 
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Initial values for the solution of equations (6) are : 

dx T /dq = (0, 0, 0, 0, 1). 

In the case of 

\B X /By\ , \B Z /By\ , \t X \ , \ty\ <^ 1 

the following relations between A x , A y and their derivatives are valid: 



(7) 



\dA XiV /dt Xt 
dA x(y) /dt x 
\A„\ 



< 1 



A x (y) 

A T 



The procedure of evaluating the derivatives is simplified when these relations are taken 
into account. 

Approximation A: The vector dx/dt x o can be approximated as : 

dx T /dt x0 = (dx/dt x0 , dy/dt x0 , 1, dt y /dt x0 , 0) , 

where the derivatives are solutions of the system of equations 

d/dz(dx/dt x0 ) = dt x /dt x0 , 

d/dz(dy/dt x0 ) = dt y /dt x0 , 

d/dz(dt x /dt x0 ) = , 

d/dz(dt y /dt x0 ) = q -v- [(dA y /dt x )(dt x /dt x0 ) + (dA y /dt y )(dt y /dt x0 )] , 

dq/dt x0 = 

with initial values (5). The solution for the derivatives dx/dt y0 in the same approximation 
is similar : 



The derivatives 



dx T /dt y0 = (dx/dty , dy/dt y0 , dt x /dt y0 , 1, 0) 



dx T /dq = (dx/dq , dy/dq , dt x /dq , dt y /dq , 1) 



are solutions of the equations (6) with the initial values given in (7). 

Approximation B: This is the most drastic simplification - the derivatives dA x ^ y /dt Xjy 
as well as A y are neglected in the corresponding equations. The system of equations (4) 
is simplified to : 



d/dz(dx/dt x0 ) - 
d/dz(dy/dt x0 ) = 
d/dz(dt x /dt x0 ) 
d/dz(dt y /dt x0 ) 
dq/dt x0 = 



dt x / dt x0 , 

Sty/ 8t X Q , 

0, 
0, 
. 



The solution of this system with initial values (5) is 

dx T /dt x0 = (s, 0, 1, 0, 0) . 
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The solution of similar equations for derivatives d/dt y0 is 



dx T /dt y0 = (0, s, 0, 1, 0) . 
For the vector of derivatives d/dq we obtain : 

dx T /dq = (dx/dqo, 0, dt x /dq , 0, 1) , 
where the derivatives are solutions of the system of equations 

d/dz(dx/dq ) = dt x /dq , 



(8) 

d/dz(dt x /dq ) = v ■ A x , 



with initial values from (7). 



4 Application of the Kalman filter technique for the 
magnet tracking 

We use the notations from [9] to describe the practical implementation of the Kalman 
filter technique for track following in the magnet. The system state vector (1) after 
inclusion of k measurements is denoted by x k , and its covariance matrix by C k . The 
coordinate measurement corresponding to the hit k is denoted by m k . The HERA-B 
magnet tracking detectors (drift tubes and Micro-Strip Gaseous Chambers) measured 
only one coordinate and m k is a scalar and its covariance matrix V k contains only one 
element. The relation between the track parameters x k and the expected measurement 
m k is described by the projection matrix H k . The matrix H k has the structure: 

= (h u 0, 0, 0, 0) , (9) 

for a detector plane measuring only x coordinates (signal wires of the drift tubes or anode 
strips of the MSGC are parallel to the y axis), and 

= (h u h 2 , 0, 0, 0), (10) 

for stereo planes rotated around the z axis. The predicted state vector x 1 ^ 1 is determined 
as the solution of the equations (2) with the initial value x^-i- The covariance matrix of 
the vector x\~ l is obtained by the propagation of the matrix C k -i'- 

C k 1 = F k C k -\F k + Q k , (11) 

where Q k denotes the covariance matrix of the process noise and the transport matrix is 

(i 2 ) 

0(a*-i) 

The estimated residual and variance become 

rt 1 = m k -H k x k k -\ 

pfc-l _ T/ , TT ^ffc-l TTT 

ti k — v k + ti k U k H k . 
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The updating of the system state vector after inclusion of the measurement k is obtained 
with the filter equations: 

K k = C k k ~ l H T {RtW 

x k = x k - 1 + K k r k '\ (14) 
C k — (1 - K k H k ) C k , 

with the filtered residuals 

r k = (l-H k K k )r k k -\ 
Rk — 0- — H k K k ) V k . 

The x 2 contribution of the filtered point is given by: 

xl = rlR k l r k . 




7.-7. . \cm\ 7-7 . \cm\ 

magnet L J magnet L J 

Figure 1: Magnetic field components as a function of the z coordinate in the magnet center (a) 
and at a vertical displacement of 40 cm (b). The z coordinate is given relative to the magnet 
center at z magnet = 450 cm. 



Because of sparse projection matrices (9,10), the calculation in (13-15) becomes rather 
simple. For the matrix the variance of estimated residuals is 

{^ _1 }n = {V k }u + h (h C u + h 2 C 21 ) + h 2 (h C 12 + h 2 C 22 ), 

i.e. it involves 6 multiplications. Here we use the notation CV,- = {C^ -1 }^-. For the same 
case, each of the five elements of the gain matrix is calculated as 

{K} a = h = (Cn h + C l2 h 2 ) {R k k ~ 1 }^ . 

The calculation of the gain matrix includes 15 multiplications and 1 division. The most 
time consuming operations are the propagation of the covariance matrix in (11) and the 
calculation of the covariance matrix in the filter equations (14) 
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A typical behavior of the magnetic field components [10] as a function of z is shown 
in fig. I. The main bending component (B y ) has a bell-shaped ^-dependence. The 
components B x , B z are sizeable away from the central axis (fig. lb). There is clearly no 
region inside the magnet where the field can be regarded as homogeneous. 

"Numerical experiments" with the real track finding procedures have shown that the 
approximation B for derivatives is accurate enough to be used for the propagation of the 
covariance matrix (11) in the case of the inhomogeneous magnetic field of HERA-B. In 
approximation B the transport matrix has a rather sparse structure: 
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(16) 



where Sk = Zk — Zk-i, and derivatives are denoted x' k = dx k k ~ x j dq^-x and t' k = dt^ 1 / dqk-i 
For a short distance (\s\ < 20 cm) the derivatives are obtained in a parabolic expansion 

as: 



dx k k 



1 / dqo = -v A x s\ and dt k xk 1 /dq = v A x s k . 



For a long distance we find the derivatives as the solution of the system (8) (initial 
values from (7)), together with the "zero trajectory" (2), by the fourth or fifth-order 
Runge-Kutta method [8]. 

The propagation of the covariance matrix F k Ck-\F k in (11) we perform in two steps. 
First, we define the product of two matrices U = F k C k -i- 



U\\ — Cn + C13 S k + C15 X k , 

Ul2 = C12 + C 23 S k + C 25 X k , 

Ul3 = C13 + C33 S k + C35 X k , 

Uu = C M + C34 S fc + C45 X k , 

Ul5 = C15 + C35 Sfc + C55 X k , 

u 21 — Cl2 + C14 Sfc, 

u 22 — C22 + C24 Sk, 

u 23 — C23 + C34 Sk, 

U24 = c 24 + C44 s k , 
U25 = C 2 5 + C45 s k , 
u 31 — c 13 + c 15 t'ki 
u 32 = C23 + C25 t' k , 
u 33 — c 33 + c 35 t'ki 
M34 = C34 + C45 t k , 
U35 = C35 + c 55 t' k , 

Uii — 



(17) 



(i 



1-5-5), 
1-5), 



where c^- = {Ck-i}ij and Uij = {U}ij. Then the matrix U is multiplied by F k to obtain 



7 



the final symmetric matrix C = F k C k ^iF k : 



C21 

C31 
C 41 

c f 

C22 
C 3 2 
C42 
C52 
C33 



c 4 

C54 

c 55 



43 
53 



41 



+ U13 S k + U 15 X k , 

U21 + u 2 3 s k + U25 x' k , 

U31 +U 33 Sfc + M35 x'k, 
U41 + M43 Sfc + M45 ^4, 
U51 + u 53 s k + m 55 x' k , 
U22 + U 2 4 s k , 
U32 + M34 s k , 
U42 + U44 s k , 
U52 + M54 Sk, 
U33+U 35 t' k , 



(18) 



U43 - 
U53 - 
U44, 

U54, 

«55- 



M 45 4) 
U55 t'k. 



The evaluation of the elements in (17) and (18) implies 37 multiplications, that is much 
smaller than the 200 multiplications needed for the case of the completely filled matrix 
F k - 

The product of the matrices K k and in (14) is 
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and for the matrix 
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For the case with the matrix the co variance matrix Ck is given by: 



i hi — 
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Only 40 multiplications are sufficient to obtain the matrix Ck in this case and 20 multi- 
plications for the option with the matrix This has to be compared with 100 mul- 
tiplications needed for the completely filled matrix H. The optimized matrix evaluation 
was implemented in functions of the ranger package related to the magnet tracking [4]. 

5 Program implementation 

The described approach for optimized integration of the equations 2 and their derivatives 
was implemented as a set of C++ functions included in the ranger package and used for 
the magnet tracking [4] and the track fit [11, 12]. All functions are of type void. In the 
following , we list the definition of common parameters of these functions in C++: 

// 

// INPUT PARAMETERS 

// 

double z_in ; // z value for input parameters 

double p_in[5] ; // vector of input track parameter (x,y,tx,ty,Q/p) 
double c_in[25] ; // covariance matrix of input parameters 
float error [2]; // desired accuracy in cm 

// error [0] for Inner Tracker region 
// error [1] for Outer Tracker region 

// 

// OUTPUT PARAMETERS 

// 

double z_out; // z value for output parameters 
double p_out [5] ; // vector of output track parameters 
double rkd[25] ; // derivatives of output parameters with respect 
// to input 

// rkd[0] deriv. of p_in[0] with respect to p_out [0] 
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// rkd[l] p_in[0] p_out[l] 
// ... 

// rkd[5] p_in[l] p.out [0] 
// ... 

// rkd[24] p_in[4] p.out [4] 

double c_out[25];// covariance matrix of output parameters 

int ierror; // error flag ( = ok, = 1 particle curls) 

The definition of the corresponding variables in FORTRAN is : 



INPUT PARAMETERS 



real*8 z_in 
real*8 p_in(5) 
real*8 c_in(5,5) 
real error (2) 



z value for input parameters 

vector of input parameters (x,y,tx,ty,Q/p) 

covariance matrix of input parameters 

desired accuracy in cm 

error (1) for Inner Tracker region 
error (2) for Outer Tracker region 



OUTPUT PARAMETERS 



real*8 z_out 
real*8 p_out(5) 
real*8 rkd(5,5) 



! z value for output parameters 
! vector of output track parameters 
! rkd(i,j) derivative of p_in(i) with respect 
c to p_out(j) 

real*8 c_out(5,5) ! covariance matrix of output parameters 
integer ierror ! error flag ( = ok , = 1 particle curls) 

A typical function call in C++ 

rk4order_(double& z_in, double* p_in,double& z_out, double* p_out, 
double* rkd, intfe ierror); 

invokes function integrating equations for particle parameters and equations for deriva- 
tives in the approximation A by a fourth-order Runge-Kutta method . 
In FORTRAN this function can be invoked like subroutine : 

call rk4order (z_in, p_in, z_out, p_out, rkd, ierror) 

The differences are evident and in following, the call statements in FORTRAN will be 
mentioned only. 
A statement 

call rk4f ast (z_in, p_in, z_out, p_out, rkd, ierror) 

effects integration of the equations for the 'zero trajectory' and calculation of derivatives 
in the approximation B by a fourth-order Runge-Kutta method . 
The function called as : 
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call rklf ast (z_in, p_in, z_out, p_out, rkd, ierror) 

calculates parameters and derivatives in approximation B using a parabolic expansion of 
the particle trajectory. 

The function evaluating particle parameters (approximation A for derivatives) by a 
fifth-order Runge-Kutta method with adaptive step size control is invoked by the state- 
ment : 

call rk5order (z_in, p_in, error, z_out, p_out, rkd, ierror) 

The corresponding function which evaluates the derivatives in approximation B is exe- 
cuted by : 

call rk5f ast (z_in, p_in, error, z_out, p_out, rkd, ierror) 

The function invoked as: 

call rk5numde(z_in, p_in, error, z_out, p_out, rkd, ierror) 

evaluates the output parameters by a fifth-order Runge-Kutta method with adaptive 
step size control and calculates derivatives by 'numerical differentiation' of the output 
parameters as a function of the input parameters. The field gradients are not neglected 
in this case but the function spends by factor of 5 more computing time than rk5f ast. 

The "fully automatic" function rktrans transports particle parameters from plane 
z to plane z by means of three different methods depending on h = z — z as it was 
described earlier. The derivatives are calculated in approximation B. The call statement 
for this function is : 

call rktrans (z_in, p_in, z_out, p_out, rkd, ierror) 

The function rktrans c uses a similar approach to transport particle parameters and 
the corresponding covariance matrix . The structure of the derivative matrix in approx- 
imation B is fully exploited for maximum speed. The function can be invoked by a 
statement 

call rktransc(z_in, p_in, c_in, z_out, p_out, c_out, ierror) 

All described functions were designed for the conditions of the complete geometry of the 
HERA-B detector where typical transport distances are 50 — 100 cm. The limitations of 
tracking precision coming from multiple scattering and measurement errors in trackers 
constrain the tracing accuracy required for such distances. 

The function rk5clip can be used to propagate particle parameters over larger dis- 
tances with higher accuracy (approximation A for derivatives). The call statement for 
this function is 

call rk5clip(z_in, p_in, z_out, p_out, rkd, ierror) 
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/im 


[im 


5 


12 


13 


10 


9 
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30 
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60 


7 


6 


90 


7 


3 



Table 1: The computational accuracy of rk5clip for different momenta 

The following procedure was used for the tuning of the steering parameters. Parti- 
cle parameters were generated with production slopes t x>y uniformly distributed from 
— 100 mrad to 100 mrad. Each particle was traced by small steps from the target to 
the area behind the magnet (z = 700 cm) and then traced back using rk5clip. The 
difference between initial and final particle coordinates was regarded as a measure of 
computational accuracy. The accuracy for selected steering parameters is shown in the 
table. As expected, the function rk5clip spends more computing time especially for low 
momentum (the dependence is roughly oc 1/po) ■ It should be noted that all described 
functions do not check if the magnetic field is defined in the region where the particle 
should be traced. The user himself must make sure that a corresponding function is 
invoked within the magnetic field and should use linear line propagation in the field-free 
case. Also, an attempt to trace a very slow particle will not be successful because the 
equations (1) do not describe a particle curling in the magnetic field. 
An additional function, called as : 

real zmin,zmax ! in cm 

call rkzf ield(zmin,zmax) 

returns as output the lower and upper bounds zmin , zmax (with respect to the center of 
the magnet) of the region where the field is defined by routines guf Id, utf eld [10]. 

This region ( roughly ±5 m from the center of the magnet) includes the area where 
we have field measurements or at least results from the MAFIA calculation. At the edges 
of the region the magnetic field can be neglected so that simple linear line searches and 
fits are sufficient for pattern recognition and track fitting [13]. This can save computing 
time for event generation and reconstruction. Therefore during the event simulation in 
HBGEAN [14] the procedure for particle tracing in the magnetic field is invoked only 
when the particle is within the GEANT volume MAGN . The definition of this volume 
can be found in the Arte table GESL for the detector component Magnet. Note that this 
volume is not identical with the region of non-zero field as defined by guf Id, utf eld. 
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